Discrete Fourier Transform and Fast Fourier Transform#
Last edited: Tom Abel Feb 2022 Tom Abel, 2021: Added painting with circles section, interpolating with FFTs. ___This Notebook gives a brief introduction to Discrete Fourier Transform (DFT) and the Fast Fourier Transform (FFT). The radix-2 Cooley-Tukey FFT algorithm is implemented and toward the end the physical meaning is explained.
These concepts have a wide area of applications in many different areas in both physics and mathematics, such as signal processing, sound and image filtering, data compression, partial differential equations and multiplication of large integers.
We explore one application in which we draw a smooth curve in the plane through a small number of input points by using the DFT to calculate interpolation coefficients.
We also will introduce spectral densities also called power spectrum sometimes. A concept that you will undobtedly run into soon if you have not already as it pops up in an enormous number of fields of study.
We start by importing the needed packages.
import numpy as np
import matplotlib.pylab as plt
from numpy import random as rnd
import timeit
from scipy import fftpack as fft
Discrete Fourier Transform (DFT)#
Let \(\vec x = [x_0,x_1,...,x_{n-1}]\) be a vector with \(n\) complex (or real) elements. The DFT of \(\vec x\) is the complex vector \(\vec y = [y_0,y_1,...,y_{n-1}]\), where the elements are defined as $\(y_k=\sum_{j=0}^{n-1}x_j\omega^{k\cdot j},\)\( where \)\omega = \exp(-2\pi i /n)\( (\)i$ is the imaginary unit) [1].
def DFT(x):
""" Calculates the one dimensional discrete Fourier transform of
a vector.
:x: double arr. The vector that is being transformed.
:returns: double arr. The Fourier transform of x.
"""
n = len(x)
y = [0]*n
omega = np.exp(-2.0j*np.pi/n)
for k in range(0,n):
y[k] = np.sum(x*omega**(np.arange(0,n)*k))
return y
It is easy to realize that the inverse DFT is given by $\(x_k = \sum_{j=0}^{n-1} y_j\omega^{k\cdot j},\)\( where \)\omega = \exp(2\pi i/n)$.
def inverseDFT(y):
""" Calculates the inverse one dimensional discrete Fourier
transform of a vector.
:y: double arr. The vector that is being transformed.
:returns: double arr. The inverse Fourier transform of y.
"""
n = len(y)
x = [0]*n
omega = np.exp(2.0j*np.pi/n)
for k in range(0,n):
x[k] = np.sum(y*omega**(np.arange(0,n)*k))/float(n)
return x
Let us try with a small example where we simply transform and inverse transform an arbitrary vector.
# Defining an array that is being transformed.
x = rnd.randint(8,size=8)
print('x =', x)
# The Fourier transform
y = DFT(x)
print('y =', np.round(y,2))
# The inverse Fourier transform.
x = inverseDFT(y)
print('x =', np.round(x,2))
x = [7 6 3 0 4 2 0 7]
y = [29. +0.j 10.78-0.88j 8. -1.j -4.78+5.12j -1. -0.j -4.78-5.12j
8. +1.j 10.78+0.88j]
x = [ 7.-0.j 6.-0.j 3.-0.j -0.+0.j 4.+0.j 2.+0.j 0.+0.j 7.+0.j]
As you already might have noticed, this DFT-algorithm is quite inefficient. There are many subcalculations that are performed more than once, and as a consequence the complexity of this algorithm is \(\mathcal O(n^2)\).
Fast Fourier Transform (FFT)#
The FFT algorithms exploits symmetries and that many operations are similar. In this notebook we are going to discuss the Cooley–Tukey algorithm [2].
Assume that \(N\) is composite. This means that \(N=n_1\cdot n_2\), where \(N\), \(n_1\) and \(n_2\) are integers. Rewrite the two indicies as $\(k=n_2k_1+k_2,\)\( \)\(j = n_1j_2 + j_1,\)\( where \)k_{1,2} = 0,1,…,n_{1,2}-1\( and \)j_{1,2} = 0,1,…,j_{1,2}-1\(. If we insert these new indicies into the DFT, some cross terms vanish, and the final result is \)\(y_{n_2k_1+k_2}=\sum_{j_1=0}^{n_1-1}\sum_{j_2=0}^{n_2-1}x_{n_1j_2+n_1}\exp\left[\frac{-2\pi i}{n_1n_2}(n_1j_2+j_1)(n_2k_1+k_2)\right]\)\( \)\(=\sum_{j_1=0}^{n_1-1}\exp\left[-\frac{2\pi i}{n}j_1k_2\right]\left(\sum_{j_2=0}^{n_2-1}x_{n_1j_2+j_1}\exp\left[-\frac{2\pi i}{n_2}j_2k_2\right]\right)\exp\left[-\frac{2\pi i}{n_1}j_1k_1\right].\)\( In this equation each inner sum is a DFT of size \)n_2\( and each outer sum is a DFT of size \)n_1\(. This yields a recursive formula for computing the DFT, which is explained in more detail in [1] and [4]. For simplicity, let us use the radix-2 algorithm. The complexity of the FFT algorithm is \)\mathcal O (n\log n)$, which makes it almost linear for large data sets!
def CooleyTukeyRadix2FFT(x):
""" Calculates the one dimensional discrete Fourier transform of
a vector using the radix-2 Cooley-Tukey FFT algorithm. The vector
that is being transformed must have a power of 2 number of elements.
:x: double arr. The vector that is being transformed.
:returns: double arr. The Fourier transform of x.
"""
# Check if n is a power of 2.
if ( len(x) & (len(x) - 1)):
raise Exception("The number of elements in x has to be a power of 2!")
# Recursive formula for calculating the FFT.
def foo(x):
n = len(x)
if n == 1:
y = x
else:
y2 = foo(x[0:n:2])
y1 = foo(x[1:n + 1:2])
d = np.exp(-2j*np.pi/n)**np.arange(0,n/2)
y = np.append(y2 + d*y1,y2 - d*y1)
return y
return foo(x)
def inverseCooleyTukeyRadix2FFT(y):
""" Calculates the one-dimensional inverse discrete Fourier transform of
a vector using the radix-2 Cooley-Tukey FFT algorithm. The vector
that is being transformed must have a power of 2 number of elements.
Parameters:
x: double arr. The vector that is being transformed.
Returns:
y: double arr. The Fourier transform of x.
"""
# Check if n is a power of 2.
if (len(y) & (len(y) - 1)):
raise Exception("The number of elements in x has to be a power of 2!")
# Recursive formula for calculating the FFT.
def foo(y):
n = len(y)
if n == 1:
x = y
else:
x2 = foo(y[0:n:2])
x1 = foo(y[1:n + 1:2])
d = np.exp(2j*np.pi/n)**np.arange(0,n/2)
x = np.append(x2 + d*x1,x2 - d*x1)
return x
return foo(y)/len(y)
Let us try with a small example where we simply transform and inverse transform an arbitrary vector as before.
# Defining the array that is being transformed.
x = rnd.randint(10,size=8)
print('x =', x)
# The Fourier transform.
y = CooleyTukeyRadix2FFT(x)
print('y =', np.round(y,2))
# The invese Fourier transform.
x = inverseCooleyTukeyRadix2FFT(y)
print('x =', np.round(x,2))
x = [8 5 2 6 9 2 5 0]
y = [37. +0.j -3.12-3.36j 10. -1.j 1.12-9.36j 11. +0.j 1.12+9.36j
10. +1.j -3.12+3.36j]
x = [8.+0.j 5.-0.j 2.+0.j 6.+0.j 9.-0.j 2.+0.j 5.-0.j 0.-0.j]
To demonstrate the superiority of the FFT we calculate the Fourier transform of a lot bigger data set. Let us also compare with the fft function from scipy.fftpack.
x = rnd.rand(8192)
# Time the loop time for DFT, CooleyTukeyRadix2FFT, numpy.fft, and scipy.fftpack.fft .
%timeit y = DFT(x)
%timeit y = CooleyTukeyRadix2FFT(x)
%timeit y = np.fft.fft(x)
%timeit y = fft.fft(x)
3.53 s ± 20.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[7], line 4
2 # Time the loop time for DFT, CooleyTukeyRadix2FFT, numpy.fft, and scipy.fftpack.fft .
3 get_ipython().run_line_magic('timeit', 'y = DFT(x)')
----> 4 get_ipython().run_line_magic('timeit', 'y = CooleyTukeyRadix2FFT(x)')
5 get_ipython().run_line_magic('timeit', 'y = np.fft.fft(x)')
6 get_ipython().run_line_magic('timeit', 'y = fft.fft(x)')
File ~/opt/anaconda3/lib/python3.9/site-packages/IPython/core/interactiveshell.py:2414, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
2412 kwargs['local_ns'] = self.get_local_scope(stack_depth)
2413 with self.builtin_trap:
-> 2414 result = fn(*args, **kwargs)
2416 # The code below prevents the output from being displayed
2417 # when using magics with decodator @output_can_be_silenced
2418 # when the last Python token in the expression is a ';'.
2419 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):
File ~/opt/anaconda3/lib/python3.9/site-packages/IPython/core/magics/execution.py:1174, in ExecutionMagics.timeit(self, line, cell, local_ns)
1171 if time_number >= 0.2:
1172 break
-> 1174 all_runs = timer.repeat(repeat, number)
1175 best = min(all_runs) / number
1176 worst = max(all_runs) / number
File ~/opt/anaconda3/lib/python3.9/timeit.py:205, in Timer.repeat(self, repeat, number)
203 r = []
204 for i in range(repeat):
--> 205 t = self.timeit(number)
206 r.append(t)
207 return r
File ~/opt/anaconda3/lib/python3.9/site-packages/IPython/core/magics/execution.py:158, in Timer.timeit(self, number)
156 gc.disable()
157 try:
--> 158 timing = self.inner(it, self.timer)
159 finally:
160 if gcold:
File <magic-timeit>:1, in inner(_it, _timer)
Cell In[5], line 23, in CooleyTukeyRadix2FFT(x)
21 y = np.append(y2 + d*y1,y2 - d*y1)
22 return y
---> 23 return foo(x)
Cell In[5], line 19, in CooleyTukeyRadix2FFT.<locals>.foo(x)
17 else:
18 y2 = foo(x[0:n:2])
---> 19 y1 = foo(x[1:n + 1:2])
20 d = np.exp(-2j*np.pi/n)**np.arange(0,n/2)
21 y = np.append(y2 + d*y1,y2 - d*y1)
Cell In[5], line 18, in CooleyTukeyRadix2FFT.<locals>.foo(x)
16 y = x
17 else:
---> 18 y2 = foo(x[0:n:2])
19 y1 = foo(x[1:n + 1:2])
20 d = np.exp(-2j*np.pi/n)**np.arange(0,n/2)
Cell In[5], line 19, in CooleyTukeyRadix2FFT.<locals>.foo(x)
17 else:
18 y2 = foo(x[0:n:2])
---> 19 y1 = foo(x[1:n + 1:2])
20 d = np.exp(-2j*np.pi/n)**np.arange(0,n/2)
21 y = np.append(y2 + d*y1,y2 - d*y1)
Cell In[5], line 19, in CooleyTukeyRadix2FFT.<locals>.foo(x)
17 else:
18 y2 = foo(x[0:n:2])
---> 19 y1 = foo(x[1:n + 1:2])
20 d = np.exp(-2j*np.pi/n)**np.arange(0,n/2)
21 y = np.append(y2 + d*y1,y2 - d*y1)
Cell In[5], line 18, in CooleyTukeyRadix2FFT.<locals>.foo(x)
16 y = x
17 else:
---> 18 y2 = foo(x[0:n:2])
19 y1 = foo(x[1:n + 1:2])
20 d = np.exp(-2j*np.pi/n)**np.arange(0,n/2)
Cell In[5], line 18, in CooleyTukeyRadix2FFT.<locals>.foo(x)
16 y = x
17 else:
---> 18 y2 = foo(x[0:n:2])
19 y1 = foo(x[1:n + 1:2])
20 d = np.exp(-2j*np.pi/n)**np.arange(0,n/2)
Cell In[5], line 21, in CooleyTukeyRadix2FFT.<locals>.foo(x)
19 y1 = foo(x[1:n + 1:2])
20 d = np.exp(-2j*np.pi/n)**np.arange(0,n/2)
---> 21 y = np.append(y2 + d*y1,y2 - d*y1)
22 return y
KeyboardInterrupt:
Physical meaning#
The DFT maps a finite equally spaced sample sequence from its original domain to its frequency domain. In other words, a discrete time data set are transformed into a discrete frequency data set.
To illustrate this, we need to figure out what the DFT-formula physically means. We start by rewriting it as $\(x_k=\sum_{j=0}^{n-1}y_j\exp\left(2\pi i\frac{k}{n\Delta t}j\Delta t\right).\)\( What the expression tells us is simply that \)\vec x\( is a superposition of exponential functions with different frequencies \)f_j = \frac{j}{n\Delta t}\( and amplitudes \)y_j\(. Therefore, we can view the magnitude of the amplitudes \)|y_k|^2\( as a measure of the "weight of the frequency \)f_j\(" in \)\vec x$!
Do remember that the complex exponential functions relate to cosines and sines through Euler’s formula
Multidimensional DFT#
Let \(\vec{j} = (j_1,j_2,...,j_d)\) and \(\vec{k} = (k_1,k_2,...,k_d)\) be \(d\)-dimensional vectors of indicies from \( 0\) to \(n-1 = (n_1-1,n_2,...,n_d-1)\). Then, the \(d\)-dimensional DFT is given by $\(y_{\vec{k}}=\sum_{\vec{j}=0}^{n-1}x_{\vec{j}}\exp\left[-2\pi\vec{k}\cdot\vec{\xi}\right],\)\( where \)\vec{\xi}\( is the elementwise division \)(j_1/n_1,…,j_d/n_d)\( [4]. For example, the two dimensional DFT is given by \)\(\vec {y_{k_1,k_2}}=\sum_{j_1=0}^{n_1-1}\sum_{j_2=0}^{n_2-1}x_{j_1,j_2}\exp\left[-2\pi i\left(\frac{ k_1j_1}{n_1}+\frac{k_2j_2}{n_2}\right)\right].\)$
Interpolation using FFTs: Drawing with circles#
This is a very fun activity and made popular by this really neat youtube video among others by 3Blue1Brown. Someone on stackexchange was working through it and drew a nice version of \(\pi\)
and stored a set of 120 x and y positions along the symbol in a handy comma separated file, which I turned into a list to include here.
xy = [108,0,110,25,112,50,113.5,75,115,100,116,125,117.5,150,125,150,150,150,175,150,200,150,225,150,225,175,225,200,225,220,200,220,175,220,150,220,125,220,100,220,75,220,50,220,25,220,0,219.5,-25,219,-50,217,-75,215,-100,212,-125,209,-150,203,-158,200,-175,190,-190,175,-203,150,-211,125,-220,100,-225,85,-209,85,-200,100,-182,125,-175,132,-150,145,-125,150,-100,150,-87,150,-87.5,125,-89,100,-92,75,-95,50,-100,25,-105,0,-113,-25,-122,-50,-136,-75,-152,-100,-170,-125,-186,-150,-189,-175,-178,-200,-175,-205,-150,-220,-125,-220,-100,-202,-85,-175,-77,-150,-73,-125,-70,-100,-67.5,-75,-65,-50,-62,-25,-60,0,-57,25,-54.5,50,-51.5,75,-49,100,-47,125,-45,150,-25,150,0,150,25,150,50,150,58,150,55,125,53,100,51,75,49,50,47,25,44.5,0,42,-25,40,-50,38.5,-75,37.5,-100,37,-125,37.5,-150,43,-175,49,-185,66,-200,75,-205,100,-215,125,-218,150,-214,175,-203,179,-200,201.5,-175,213,-150,221,-125,226.5,-100,227.5,-88,210,-88,209,-100,200,-123,197,-125,175,-141,150,-144,125,-134,117,-125,109,-100,106,-75,106,-50,106.5,-25]
x = np.array(xy[::10]) # select only a subset of the values
y = np.array(xy[1::10])
N = len(x)
plt.plot(x,y);
So we have a parametric curve where a vectors x and y coordinate are given along the line parameterized by an extra parameter \(t\), specifying \(x(t)\) and \(y(t)\).
from scipy.interpolate import interp1d
t = np.arange(N)/N
xt = interp1d(t, x, bounds_error=False)
yt = interp1d(t, y, bounds_error=False)
xtc = interp1d(t, x, kind="cubic", bounds_error=False)
ytc = interp1d(t, y, kind="cubic", bounds_error=False)
tf = np.linspace(0,1,1200) # finer representation of t
plt.plot(t, x,'.', label="x")
plt.plot(t, y,'.', label="y")
plt.plot(tf, xtc(tf), label="x interp")
plt.plot(tf, ytc(tf), label="y interp");
plt.legend();
plt.plot(xtc(tf), ytc(tf))
[<matplotlib.lines.Line2D at 0x132460df0>]
We can interpret every point as a location in the complex plane so that any of our input \((x_j,y_j)\) points can be expressed as a complex number \(z_j = x_j + i y_j\). The inverse Fourier transform of these complex numbers will describe properties of cyclic functions which when all added together will give back the input values. So we can get using our discrete Fourier transform the values $\( Z_k = \sum_{n=0}^{N-1} z_j e^{-i jk\frac{2\pi}{N}} \)\( and can get back the values \)z[j]\( if we have \)Z[k]\( with the inverse of that transform \)\( z_j = \frac{1}{N}\sum_{k=0}^{N-1} Z_k e^{i jk\frac{2\pi}{N}} \)\( For the discrete Fourier transform most implementations consider though the form in which both positive and negative frequencies are considered: \)\( Z_k = \sum_{n=-N/2+1}^{N/2} z_j e^{-i jk\frac{2\pi}{N}} \)\( and can get back the values \)z[j]\( if we have \)Z[k]\( with the inverse of that transform \)\( z_j = \frac{1}{N}\sum_{k=-N/2+1}^{N/2} Z_k e^{i jk\frac{2\pi}{N}} \)\( Now if we think of our distinct positions \)j\( as part of a continuous field \)j \equiv t N\( then we can represent \)\( z(t) = \frac{1}{N}\sum_{k=-N/2+1}^{N/2} Z_k e^{i 2\pi k t} \)\( which will match \)z_j\( when \)t=j/N\(. This is the key idea behind Fourier interpolation. Below we will use \)z_k\equiv Z_k/N\( so that \)z(t)=\sum z_k e^{i 2\pi k t} $.
z = x[:] + 1j * y[:] # express our points on the curve as complex numbers
N = len(z)
Zk = np.fft.fft(z) # calculate Fourier coefficients
zk = Zk/N # absorb the 1/N factor into the Zk to
k = np.fft.fftfreq(N, d=1/N); # use helper function to tell us where which frequencies are
amp = np.abs(zk)
phase = np.angle(zk)
tp = 0.78
testp = np.sum(zk*np.exp(1j*k*2*np.pi * tp))
fig = plt.figure(figsize=(10,10))
ax = plt.gca()
ax.cla() # clear things for fresh plot
#ax.plot(xt(tf),yt(tf))
ax.plot(xtc(tf),ytc(tf))
ax.plot([np.real(testp)], [np.imag(testp)],'X')
xp = 0
yp = 0
xl = []
yl = []
for i in range(N):
xl.append(xp)
yl.append(yp)
# print(k[i],amp[i],phase[i], xp, yp)
# xp = xp + amp[i]*np.cos(2.*np.pi*k[i]*tp + phase[i] )
# yp = yp + amp[i]*np.sin(2.*np.pi*k[i]*tp + phase[i] )
zc = zk[i]*np.exp(1j*k[i]*2*np.pi * tp) # this is the same as above
circ = plt.Circle((xp, yp), amp[i], color='black', fill=False, alpha=.3)
plt.arrow(xp,yp,zc.real,zc.imag,overhang=.8,color="red",\
head_width=4,length_includes_head=True,alpha=.5)
xp = xp + zc.real
yp = yp + zc.imag
ax.add_patch(circ)
circ = plt.Circle((xp, yp), amp[i], color='black', fill=False, alpha=.3)
ax.add_patch(circ)
xl.append(xp)
yl.append(yp)
#ax.plot(xl,yl, color="red",alpha=.3);
Let’s try and also do an animation of a similar plot.
from matplotlib import animation
from IPython.display import HTML
def Circles2DMovie(ik,izk):
Nt = 240 # frames
print("Nt:",Nt)
amp = np.abs(izk)
phase = np.angle(izk)
k = ik
zk = izk
xyt = np.zeros((Nt,2))
for i in range(Nt):
zc = (zk * np.exp(1j*k*2*np.pi * i/Nt)).sum()
xyt[i,0] = zc.real
xyt[i,1] = zc.imag
# create a simple animation
fig = plt.figure(figsize=(10,10))
ax = plt.axes(xlim=(-250, 250), ylim=(-250, 250))
line, = ax.plot([], [], 'o',lw=5)
uline, = ax.plot([], [], '-',lw=3.3,alpha=.4)
qax = ax.quiver(zk.real,zk.imag,zk.real,zk.imag,alpha=.3)
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.close() # this makes sure the empty initial plot is not shown separately
def init():
line.set_data([], [])
return line, #line2,
def animate(i):
# ax.clear()
zp = (zk*np.exp(1j*k*2*np.pi * i/Nt))
zc = zp.cumsum()
zo = np.roll(zp.cumsum(),1)
zo[0] = 0+0j
qax.set_offsets(np.c_[zo.real,zo.imag])
qax.set_UVC(zp.real,zp.imag)
# would be nice to animate the circles [but need to figure out how to update them ..]
# for ii,czp in enumerate(zc):
# ax.add_artist(plt.Circle((czp.real,czp.imag), radius=np.abs(zp[ii]), color='black', alpha=0.3, fill=False))
line.set_data(zc[-1].real,zc[-1].imag)
uline.set_data(xyt[0:i,0],xyt[0:i,1])
return line, uline, #line2,
return animation.FuncAnimation(fig, animate, init_func=init,
frames=Nt, interval=60, blit=True, repeat=True)
mov = Circles2DMovie(k,zk)
HTML(mov.to_jshtml())
Nt: 240
/var/folders/8h/csrqy1dd1x1d93qvzyh33kn40000gn/T/ipykernel_45710/1637310528.py:49: MatplotlibDeprecationWarning: Setting data with a non sequence type is deprecated since 3.7 and will be remove two minor releases later
line.set_data(zc[-1].real,zc[-1].imag)